Nigelladine A among Selected Compounds from Nigella sativa Exhibits Propitious Interaction with Omicron Variant of SARS-CoV-2: An In Silico Study

COVID-19 has been a threat to the entire world for more than two years since its outbreak in December 2019 in Wuhan city of China. SARS-CoV-2, the causative agent, had been reported to mutate over time exposing new variants. To date, no impeccable cure for the disease has been unveiled. This study outlines an extensive in silico approach to scrutinize certain phytochemical compounds of Nigella sativa (mainly the black cumin seeds) targeting the spike protein and the main protease (Mpro) enzyme of the Omicron variant of SARS-CoV-2. The objective of this study is to investigate the extracted compounds with a view to developing a potential inhibitor against the concerned SARS-CoV-2 variant. The investigation contemplates drug-likeness analysis, molecular docking study, ADME and toxicity prediction, and molecular dynamics simulation which have been executed to elucidate different phytochemical and pharmacological properties of the tested compounds. Based on drug-likeness parameters, a total of 96 phytochemical compounds from N. sativa have been screened in the study. Interestingly, Nigelladine A among the compounds exhibited the highest docking score with both the targets with the same binding affinity which is −7.8 kcal/mol. However, dithymoquinone, kaempferol, Nigelladine B, Nigellidine, and Nigellidine sulphate showed mentionable docking scores. Molecular dynamics up to 100 nanoseconds were simulated under GROMOS96 43a1 force field for the protein-ligand complexes exhibiting the top-docking score. The root mean square deviations (RMSD), root mean square fluctuations (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), and the number of hydrogen bonds have been evaluated during the simulation. From the findings, the present study suggests that Nigelladine A showed the most promising results among the selected molecules. This framework, however, interprets only a group of computational analyses on selected phytochemicals. Further investigations are required to validate the compound as a promising drug against the selected variant of SARS-CoV-2.


Introduction
has been a trending issue since its outbreak in late December 2019 in the city of Wuhan, Hubei Province, China. Te World Health Organization (WHO) announced this disease as a global pandemic on March 11, 2020 [1]. Due to the hasty spread of the virus, the socioeconomic condition of the entire world started to collapse. Te World Health Organization (WHO) reported more than 617 million cases worldwide until September 30, 2022 [2]. In the course of time, more than 6.53 million deaths occurred due to COVID-19, which is massive from the usual perspective [2]. To date, an impeccable cure for this disease is still to be unveiled.
Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), a ribonucleic acid (RNA) virus, is the causative agent of COVID-19 [3]. Undergoing extensive mutations, the virus generated a number of new variants, namely Omicron, Alpha, Beta, Delta, Gamma, and more [4]. WHO classifed these variants into two types: variants of concern (VOCs) and variants of interest (VOIs) [4]. According to WHO, the Omicron B.1.1.529 variant was frst reported to WHO on November 24, 2021, and WHO categorized it as a VOC for the frst time on November 26, 2021 [4]. Since then, this variant has been ranked above other VOCs. Tis makes Omicron B.1.1.529 a more appropriate issue to scrutinize. Tus, this study aimed to investigate if there is any promising cure that can combat this variant.
Current fndings demonstrate that mutations of SARS-CoV-2 variants are found to be more prevalent in the spike protein of the virus [5]. Hence, we have selected the spike protein as one of the targets in this study. Tis protein remains as a trimer on the surface of the viral envelope [6]. Te receptor-binding domain (RBD) is possessed by the S1 domain and is particularly liable for binding the virus to the receptor [7,8]. On the other hand, HR1 and HR2 are contained in the S2 domain, which is afliated with viral fusion [8]. Te RBD of the spike protein interacts with the host cells, admitting the receptor angiotensin-converting enzyme 2 (ACE2) [9]. Here, ACE2 itself also works as a receptor for the SARS-CoV-2 spike protein. By binding with ACE2, the virus facilitates endosome formation [9]. Consequently, it triggers viral fusion at a lower pH value [9]. Tis circumstance implies that, by intervening in the interaction between the spike protein and its receptor, the activity of the virus can be inhibited.
Te main protease (M pro ) in all variants of SARS-CoV-2 is an enzyme of SARS-CoV-2 that is broadly targeted by researchers. Tis enzyme plays an essential role in infuencing viral replication and transcription [10]. Te key function of this enzyme is to release functional polypeptides from each polyprotein through vast proteolytic processing [11]. Te functions stimulate the replication of the virus, which is the key factor in its recurrence. M pro is a homodimer that comprises two protomers each and incorporates three domains, namely domains I, II, and III [12]. Moreover, the human body does not possess a protein or enzyme nearly homologous to M pro . Tese cases note this enzyme as an ideal drug target to study. Interceding the process of this enzyme might lead to a solution to limit the recurrence of the virus. Tus, M pro is included as one of the targets in this study.
Traditional medicines are playing a pivotal role in treating miscellaneous diseases, including a notable number of viral ones [13]. According to WHO, roughly 80% of the world's population relies on traditional medicines [14]. Nigella sativa (the black cumin seeds), belonging to the family Ranunculaceae, is one of the noteworthy plants with a signifcant medicinal profle [15]. Current literature states that N. sativa was confrmed to show antiasthmatic, anticancer, anti-infammatory, antimicrobial, antioxidant, bronchodilator, hepato-protective, immunomodulator, renal protective, and many remarkable properties [16,17]. More importantly, N. sativa showed promising activities against the SARS-CoV-2 Wuhan variant, which makes this plant a noteworthy herb to inspect [17]. Recent studies explored Nigellidine, nigellicine, nigellimine, thymol, α-hederin, thymoquinone, dithymoquinone, hederagenin, etc. compounds from this plant to inhibit selective targets of coronavirus [18]. Interestingly, the activities of this plant against Omicron variant were unexplored and missing in the current literature to date. Tus, N. sativa was picked in this study to scrutinize its activity against the targets of Omicron variant SARS-CoV-2.
Computer-aided drug design (CADD) via in silico methods is a convenient approach that accelerates the process of drug discovery and development [19]. In this computational study, we employed molecular docking, pharmacokinetic and pharmacodynamic property analyses, and molecular dynamics simulations to fnd the most suitable drug candidate. Molecular docking equipped us with data about the binding afnity, orientation, and type of interactions of each ligand with the respective target proteins. Te pharmacokinetic profles were acquired to investigate the data on the absorption, distribution, metabolism, and excretion (ADME) of the compounds that occur inside the body after drug administration. Te toxicity study was carried out to get the LD50 values and toxicity classes of the individual ligands. Finally, molecular dynamics simulations were conducted to determine the stability and fexibility, along with certain properties, of the proteinligand complexes.
In the present study, we approached a screening and 96 phytochemical compounds from N. sativa were sorted out based on the drug-likeness parameters. Trough computational analyses, the molecules were tested in a wide spectrum. Hence, this framework interprets and presents some promising drug candidates from N. sativa against the SARS-CoV-2 Omicron B.1.1.529 variant.

Selection and Preparation of Ligands.
On the basis of drug-likeness, a total of 96 phytochemical molecules from Nigella sativa were selected for this study. Lipinski's rule of fve and Ghose's rules were considered while selecting the molecules [20,21]. Only the molecules following both rules were picked for the study. Te 3-dimensional (3D) conformers of the selected ligands were downloaded in SDF formats from the online databases of PubChem (https:// pubchem.ncbi.nlm.nih.gov/) and IMPPAT 2.0 (Indian Medicinal Plants, Phytochemistry, and Terapeutics; https:// cb.imsc.res.in/imppat/) [22,23].

Retrieval and Preparation of Target Protein.
Te crystal structures of the spike protein (PDB ID: 7QNW) and M pro (the main protease; PDB ID: 7TVX) of the SARS-CoV-2 Omicron B.1.1.529 variant were downloaded in PDB format from the database of the RCSB Protein Data Bank (https:// www.rcsb.org/) [24,25]. Te resolutions of the downloaded spike protein and the M pro were 2.40Å and 2.094Å, respectively. Te protein structures were cleaned by removing undesired atoms and molecules (including ligands) using PyMOL version 2.5.2 software (Schrödinger, LLC) [26]. Te receptor-binding domain (RBD) of the spike protein was isolated from the crystal structure, and the excessive chains of proteins were removed using PyMOL. Tis method was employed on M pro also, and only one protein chain was kept. Te energies of the selected protein chains of the spike protein and M pro were minimized in Swiss-PdbViewer version 4.1.0 software using preset parameters [27]. Te chains of the minimized proteins were saved in PDB formats for molecular docking and molecular dynamics simulations.

Molecular Docking.
Molecular dockings on the selected ligands were performed against the target proteins using the CB-Dock2 server (https://cadd.labshare.cn/cb-dock2/php) [28]. Te binding afnity (kcal/mol) for each protein-ligand complex as well as the noncovalent interactions and docking orientations were scrutinized by visualizing in the BIOVIA Discovery Studio 2021 Client version 21.1.0 software (Dassault Systèmes). Te schematic illustrations of the protein-ligand docking complexes were retrieved in 2D and 3D forms from BIOVIA Discovery Studio.

ADME and Toxicity Prediction.
Te canonical SMILES of the ligands with the top-docking scores were copied from the PubChem and IMPPAT 2.0 databases and were inputted on the SwissADME server (https://www.swissadme.ch/) [29]. Te ADME (absorption, distribution, metabolism, and excretion) data for each ligand were obtained from Swis-sADME. Subsequently, the toxicity profle of each ligand was predicted from the ProTox-II server (https://tox-new. charite.de/protox_II/) [30]. Te physicochemical, pharmacokinetic, and pharmacodynamic properties of each ligand were noted from these two sources. Te topological polar surface area (TPSA), lipophilicity (MLogP), water solubility (LogS), bioavailability score, blood-brain barrier (BBB) permeability, interaction with P-glycoprotein (P-gp), LD50 value, and toxicity class of each ligand were investigated during the ADME and toxicity prediction.

Molecular Dynamics Simulation. Molecular dynamics
(MD) provides data on the stability and fexibility of proteinligand complexes. Te molecular dynamics were simulated using Groningen Machine for Chemical Simulations (GROMACS) software [31]. All simulation processes were carried out using the GROMOS96 43a1 force feld. Te MD simulation was conducted for up to 100 nanoseconds for each protein-ligand complex. At the onset, ligand topology fles for each ligand were generated from the PRODRG server (https://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg/) by separating the ligands from the docked complexes retaining the same conformations and orientations [32]. During the simulation, the box was solvated with SPC water models, and the box type was set to a triclinic shape. To neutralize the system, 0.15 M NaCl salt was added to it. Structural optimization of 5000 steps was done by minimizing the energy of the system using the steepest descent algorithm. Te simulation was carried out with equilibriumtype NVT and NPT. Te temperature and the pressure of the system were maintained at 310 K and 1.0 bar, respectively, during the processes of simulation. Te MD integration was done using the leap-frog method. From the results, the root mean square deviations (RMSD), root mean square fuctuations (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), and the number of hydrogen bonds in the protein-ligand complexes were scrutinized to get the most suitable molecule.

Molecular Docking and Noncovalent Interactions
Analysis. Nigelladine A showed the best docking scores against both target proteins, with a binding afnity of −7.8 kcal/mol in both cases. Interestingly, the top 6 molecules on the basis of binding afnity remain the same for both target proteins, only varying in binding afnities. Against the spike protein, the other molecules showing good binding afnity are kaempferol, Nigellidine, dithymoquinone, Nigellidine sulphate, and Nigelladine B, of which the binding afnities are −7.6, −7.5, −7.5, −7.4, and −7.3 kcal/mol, respectively. Te binding afnities and the noncovalent interactions (hydrogen bonds and hydrophobic) of the top 6 molecules are presented in Table 1. Te 3D and 2D docked conformations of the ligands with M pro are shown in Figure 1.
Against the spike protein, Nigelladine A and Nigellidine sulphate each formed only one hydrogen bond with the residue LEU368 and showed four hydrophobic interactions with residues PHE342, PHE374, PHE375, and TRP436. Nigelladine B formed bonds with the same residues as Nigelladine A and Nigellidine sulphate had a hydrophobic interaction with LEU371. Kaempferol formed fve hydrogen bonds and exhibited one hydrophobic interaction (Table 1). Nigellidine and dithymoquinone both did not form any hydrogen bonds. Tey showed seven and two hydrophobic interactions, respectively (Table 1).
Against the main protease (M pro ), the molecules alongside Nigelladine A exhibiting good afnities are Nigellidine sulphate, Nigellidine, kaempferol, dithymoquinone, and Nigelladine B with binding afnities of  Figure 2. Nigelladine A showed no hydrogen bonds but hydrophobic bonds with residues PRO293 and PHE294. Nigellidine, Nigellidine sulphate, and dithymoquinone showed six, four, and two interactions, respectively, with half of the interactions being hydrogen bonds in each set of interactions. Nigelladine B showed two interactions, both of which are hydrophobic, and kaempferol showed three hydrogen bonds and one hydrophobic. Te interacting residues are mentioned in Table 1.

ADME and Toxicity
Analyses. Te physicochemical, pharmacokinetic, and pharmacodynamic properties of the drug have been scrutinized using the data retrieved from SwissADME and Protox-II. ADME data of the molecules show detailed information, including molecular weight (MW), topological polar surface area (TPSA), lipophilicity (MLogP), water solubility (LogS), gastrointestinal absorption, bioavailability score, blood-brain barrier (BBB) permeability, and a complete profle of the molecules. Tese parameters indicate how suitably the molecules (entering the human body) will be absorbed, distributed, metabolized, and fnally excreted. Moreover, the toxicity level of the molecules was taken into account, since no toxic molecule would be a suitable cure. Te top 6 molecules based on molecular docking scores were considered for ADME and toxicity analyses. All molecules showed zero violations of Lipinski's rule of fve and Ghose's rule. Molecules with <90Å 2 TPSA tend to be permeable to the BBB, whereas those with >140Å 2 TPSA tend to be poorly permeable to the cell membrane [33,34]. Nigelladine A and Nigelladine B exhibit the minimum TPSA, which is 29.43Å 2 , whereas kaempferol occupied the highest value with an area of 111.13Å 2 . Among the rest, three molecules showed TPSA of less than 90Å 2 ( Table 2). Tis indicates that except for kaempferol and Nigellidine sulphate, other molecules were found to be BBB permeable. All the molecules having TPSAs <140Å 2 are able to permeate the cell membrane.
Te gastrointestinal absorption of all molecules was high. Alongside, all molecules occupied the same bioavailability score, which is 0.55. However, only Nigellidine and Nigellidine sulphate were P-glycoprotein (P-gp) substrates. In terms of toxicity analysis, kaempferol and dithymoquinone took a safer place with lethal dose 50 (LD50) values of 3919 and 2300 mg/kg, respectively, belonging to toxicity class 5 (2000 < LD50 ≤ 5000) ( Table 2). Nigelladine A and Nigelladine B showed the same and the lowest value, which is 900 mg/kg, addressing toxicity class 4 (300v< LD50 ≤ 2000).
To be administered orally, the drug candidate must follow certain criteria. Hence, the bioavailability radar of Nigelladine A depicts an overview of some major physicochemical criteria (Figure 3). Lipophilicity (in terms of XLogP3) between −0.7 and 5.0, size between 150 and 500 g/ mol, TPSA between 20 and 130Å 2 , LogS (ESOL) between −6 and 0, insaturation (fraction Csp3) between 0.25 and 1, and number of rotatable bonds between 0 and 9 are favorable for proper oral bioavailability. Nigelladine A remains in the favorable zone, occupying the suitable ranges for these criteria.

Molecular Dynamics Simulation Studies.
In molecular dynamics simulation, the root mean square deviations (RMSD), root mean square fuctuations (RMSF), radius of gyration (Rg), solvent accessible surface area (SASA), and the number of hydrogen bonds in each protein-ligand complex were scrutinized to determine the most suitable drug candidate from the selected ligands.

Root Mean Square Deviations (RMSD).
Te root mean square deviations (RMSD) were evaluated to understand the stability of the protein-ligand complexes. A lower RMSD value is always favorable since it does not have any absolute rule. Te average RMSD value for ligand Nigelladine A in the complex with spike protein was ∼5.6Å ranging between ∼2.5-11.8Å, and that of the backbone of spike protein was ∼4.0Å ranging between ∼2.5-5.0Å (Figure 4). Te values for the ligand deviated unduly during the frst half of the simulation. Te latter half showed a smaller deviation. Te average RMSD of Nigelladine A in the complex of the molecule with M pro was ∼6.2Å ranging between ∼1.2-12.6Å ( Figure 5). Te outlying values in the graph peak near 95 ns in the simulation. Te RMSD value of the backbone in the complex of Nigelladine A with M pro ranges between ∼1.4-4.1Å with an average of ∼3.0Å. Nigelladine A depicts a more favorable graphical plot with M pro than that with the spike protein.

Root Mean Square Fluctuations (RMSF).
Te root mean square fuctuations (RMSF) of the protein-ligand complexes were scrutinized to elucidate the fexibility of the protein structure (Figures 4 and 5). An RMSF value of 3.4Å or below is considered ideal [35]. Te RMSF value for the spike protein in the complex with Nigelladine A was observed between ∼0.93-9.45Å with an average of 2.25Å. In contrast, the residues of M pro fuctuated between ∼0.72-6.28Å with an average of 1.79Å. Hence, both the target proteins are found to fuctuate within the ideal RMSF value. Hence, M pro was observed to fuctuate less in comparison to the spike protein. Te residues of the spike protein fuctuated below 6.00Å except for residues 333, 334, and 477 ( Figure 4). Moreover, all residues of M pro except 1, 169, 276, and 300 fuctuated below 3.80Å ( Figure 5).    (Figures 4 and 5). Tere was a maximum of  a single hydrogen bond with the spike protein recorded during the frst half of the simulation. Several maxima were depicted in the graphical plots during the fnal half-time. In contrast, during the simulation with M pro , the maximal value was soon reached and repeated with diferent intervals of time.

Conclusion
Tis study aimed to investigate the promising molecules from N. sativa with a view to fnding potential inhibitors against the Omicron variant of SARS-CoV-2. Among the tested compounds, Nigelladine A demonstrated the most promising results against both target proteins. In terms of binding afnity, Nigelladine A exhibited the top scores in both cases while docking with the two targets. Likewise, in the molecular dynamics simulation, this molecule retained pertinent results. However, this study only elucidates the in silico properties and profles of the selected phytochemical compounds from N. sativa. Further experimental validation is required to confrm the activity of Nigelladine A as a potential inhibitor against the SARS-CoV-2 Omicron variant as well as for other variants. Hence, this study proposes Nigelladine A as a promising drug candidate showing favorable interactions against the studied targets of SARS-CoV-2.

Data Availability
All data are available within the manuscript.